The three-dimensional structure of a protein determines its biological function (see Berg J M, Tymoczko J L, Shyer L., Biochemistry, 5th edition, W H Freeman, 2002). Therefore, determining said structure is of vital importance for designing drugs, among others.
One way of determining the mentioned structure is by means of X-ray crystallization or nuclear magnetic resonance techniques (see P. Narayanan, Essentials of Biophysics, Second Edition, Anshan Ltd., 2010). Nevertheless, these techniques are not without imperfections and do not allow easily performing the opposite process, i.e., determining the amino acid chain which allows obtaining a certain structure (see Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006). These facts, along with increased computer processing capacity, have favored the development of molecular dynamics which consists of simulating the protein folding process or the interaction process of proteins with other molecules using a computer in recent years.
It has been demonstrated (see P. L. Freddolino, K. Schulten, Common Structural Transitions in Explicit-Solvent Simulations of Villin Headpiece Folding, Biophysical Journal Volume 97 October 2009 2338-2347, P. L. Freddolino, F. Liu, M. Gruebele, K. Schulten, Ten-Microsecond Molecular Dynamics Simulation of a Fast-Folding WW Domain, Biophysical Journal, Volume 94, Issue 10, L75-L77, 15 May 2008, S. K. Sadiq, G. De Fabritiis, Explicit solvent dynamics and energetics of HIV-1 protease flap-opening and closing, Proteins: Structure, Function, and Bioinformatics Volume 78, Issue 14, pages 2873-2885, 1 Nov. 2010) that approximation by means of Newton's mechanics is valid for modeling the dynamics (folding process or interaction process with other molecules) of large molecules such as proteins, which is known as molecular dynamics. Thus, by taking into account Newton's second law and that, in that case each atom is considered as a body of mass m, the equation determining the acceleration {right arrow over (a)} to which the latter is subjected due to the existence of forces {right arrow over (F)}i between the latter and the other atoms of the actual protein or of the solvent, is the following:
                                          ∑                          i              =              1                        N                    ⁢                                    F              →                        i                          =                  m          ·                      a            →                                              (        1        )            N being the number of atoms.
The acceleration {right arrow over (a)} to which a body is subjected is the derivative of the body velocity {right arrow over (v)} over time. Therefore, said velocity can be obtained from integrating the mentioned acceleration {right arrow over (a)} over time. Likewise, the position {right arrow over (x)} is obtained as the integral of the velocity {right arrow over (v)} over time given that the velocity is the derivative of position. Therefore, the position of each atom can be obtained by means of integrating from velocity and the latter is in turn obtained by means of integrating the acceleration which, with the forces acting on said atom known, provides the preceding equation 1.
Likewise, to simulate the actual conditions of proteins and of biomolecules in general, the molecule is surrounded with a solvent, water, simulating cell conditions. This solvent is necessary for certain biological processes such as protein folding, where the hydrophilic or hydrophobic nature is key for determining its structure.
This requires simulating, in addition to the protein under study, thousands of water molecules simulating the solvent, so the systems under study have thousands of atoms to be simulated.
For similar reasons, the so-called periodic boundary conditions are introduced to simulate experimental conditions of macroscopic systems (see T. Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Interdisciplinary Applied Mathematics series, vol. 21. Springer: New York, pp. 272-6. 2002). The system (the most solvent biomolecules) to be simulated is surrounded with periodic copies of itself in 3 directions.
This is a computationally approachable way for simulating the conditions of systems having a high number of molecules (of the order of Avogadro's number).
An efficient way of calculating the force acting on an atom due to its interaction with these copies is the so-called Particle Mesh Ewald (T. Darden, D. M. York, L. G. Pedersen, Particle mesh Ewald: An N log(N) method for Ewald sums in large Systems, J. Chem. Phys. (Communication) 98, 10089-10092 (1993)).
Before going into detail about the different types of force and the problem they entail, it is convenient to mention the basis of molecular dynamics, the solution of Newton's Second Law, equation 1.
It is a second-order differential equation with respect to the position and is therefore solved by integration.
This equation is equivalent to a first-order system of equations called Hamilton's equations, in what is called Hamiltonian formalism (see B. Leimkuhler, S. Reich, Simulating Hamiltonian Dynamics, Cambridge University Press 2004). Given that these solutions cannot be solved analytically in cases such as molecular dynamics, a numerical solution of said differential equations is necessary.
Said solution, called numerical integration, is based on going from a continuous time to a discrete time. The gap between each of the discrete time instants is called a timestep. Starting from initial position and velocity conditions, the values thereof in the entire subsequent instants are to be found. The function which allows obtaining a value for the following instant from the position and velocity value in one instant is called an integrator. Therefore, the integrator and the timestep determine simulation quality.
Given that an approximation is performed when going from a continuous time to a discrete time, the greater the timestep, the worse the solution of the differential equation.
If an excessively large timestep is chosen, the solution obtained may have nothing to do with the actual solution and it may diverge from it. In that case, the integrator is said to be unstable for that timestep. The range of timestep values for which an integrator is stable is called a stability zone (see E. Hairer, S. P. Nørsett, Gerhard Wanner, Solving ordinary differential equations I: Nonstiff problems, second edition, Springer Verlag, Berlin, 1993). The forces {right arrow over (F)}i to which each atom is subjected are of a different nature in terms of their magnitude, their variability over time and are the function of the position of the remaining atoms of the protein or of the atoms of the solvent.
The aforementioned numerical integration process is based on gradually moving the simulation forward at time intervals until reaching the end time thereof. Due to the very high frequency (variability over time) of some of the forces {right arrow over (F)}i these time intervals are very small, of the order of femtoseconds (10−15 seconds). The total simulation times required for simulating a complete protein folding, for example, range from microseconds to several minutes (David Sheehan, Physical Biochemistry, Second Edition, John Wiley & Sons, 2009), therefore the number of time intervals necessary for completing the simulation range from 1010 to 1016. In each of these time intervals, it is necessary to perform complex computational calculations, mainly those corresponding to the forces {right arrow over (F)}i, and furthermore, these calculations increase exponentially with the complexity of the system. This means that even the most powerful computers today need a very long computing time of the order of days, months and even years to complete the simulation of medium-sized molecular systems (see Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006. J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chypot, R. D. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD, Journal of Computational Chemistry, 26:1781-1802, 2005).
It is easily deduced from this that the techniques known in the state of the art in molecular dynamics have the technical problem that they do not allow performing simulations with the speed demanded by pharmaceutical industry or research organizations.
Among the mentioned {right arrow over (F)}i, the binding forces model very high frequency vibrations and therefore really need very small time intervals for being integrated: from 1 to 10 femtoseconds. In contrast, among these forces {right arrow over (F)}i there are other forces such as Van Der Waals or electrostatic forces (Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006) which, compared with the previous binding forces, are long-range forces having a much smaller frequency, and for which, in principle, much greater time intervals (of the order of two orders of magnitude more than those corresponding to binding forces) may be sufficient for being integrated. It must be highlighted that compared with binding forces which only take into account the interactions between one atom and its closest neighbor, these long-range forces have a very high computing cost since they involve the interactions of one atom with the other atoms of the protein or molecule, as well as with the atoms of the solvent. Therefore, while the evaluation of such forces has a complexity of O(N), for some the complexity is O(N2) (Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006).
Due to its simplicity and exceptional stability in long-term simulations such as those required in molecular dynamics, the integrator is the Leapfrog/Verlet/Stormer method (see L. Verlet, Computer ‘experiment’ on classical fluids: I. Thermodynamical properties of Lennard-Jones molecules, Phys. Rev., 159 (1):98-103, July 1967), commonly known in short as the Verlet method. It is an explicit, second-order integrator with the particularity that it only requires evaluating the forces once every iteration, whereas other integrators of the same order (e.g., second-order Runge Kutta) require 2 evaluations, and therefore virtually twice the computing cost.
Like the differential equations that it attempts to integrate, Verlet, as a numerical integrator, has important mathematical properties such as symplecticity and time symmetry. Therefore it more faithfully represents the system that it is simulating, maintaining certain physical properties such as energy conservation and time reversibility, and mathematical properties such as physical volume conservation, which will result in a greater stability. There are variants of this method such as the velocity Verlet or the position Verlet which are widely used and which also have these important properties.
Using symplectic integrators is common in various fields of mechanics (see B. Gladman, M. Duncan, J. Candy, Symplectic integrators for long-term integrations in celestial mechanics, Celest. Mech. Dynam. Astron., 52, 221-240 1991) due to the discovery of their greater efficiency with respect to non-symplectic integrators (see M. P. Calvo, J. M. Sanz-Serna, The development of variable-step symplectic integrators, with application to the two-body problem, SIAM J. Sci. Comput., 14, 936-952. 1993). The use of a variable timestep for a symplectic integrator is also known to reduce its efficiency (R. D. Skeel, C. W. Gear, Does variable step size ruin a symplectic integrator? Physica D, 60, 311-313. 1992) and its stability (see R. D. Skeel, Variable step size destabilizes the Störmer/leapfrog/Verlet method, BIT, 33, 172-175. 1993).
Due to the complexity and therefore the chaotic behavior of the systems simulated in molecular dynamics, with thousands of degrees of freedom, the tracking of a particular atom is not relevant. One characteristic of chaotic systems is their sensitivity to initial conditions, such that a small variation in one of the dynamic variables leads to a completely different trajectory. In a molecular system there is always an initial inaccuracy in the initial positions and velocities of the system due to temperature and to problems of measuring the positions of the atoms, so a particular trajectory will be not representative. However, statistical variables as well as the configurations corresponding to local minima of the energy, which will be the conformations that the trajectories are inclined to have, are of interest. In the case of a protein, it will correspond to its native state.
Therefore, in molecular dynamics high precision in the trajectories is not necessary, which allows using non-symplectic integrators and using a variable timestep for improving efficiency (see S. J. Stuart, J. M. Hicks, M. T. Mury, An Iterative Variable-timestep Algorithm for Molecular Dynamics Simulations, Molecular Simulation, Volume 29, Issue 3 2003, pages 177-186).
Both the aforementioned techniques and those which will be described below of imposing restrictions on the high-frequency components of binding forces, and the techniques based on increasing the interval of integration of slow forces, are modifications or variations of the mentioned Leapfrog/Verlet/Stormer method.
Among the techniques of imposing restrictions (constraints) for damping high frequencies include those known as SHAKE (see J. P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, Numerical integration of the Cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes, J. Comput. Phys., 23:327-341, 1977). RATTLE (see H. C. Andersen. Rattle: a ‘velocity’ version of the SHAKE algorithm for molecular dynamics calculations, J. Comput. Phys., 52:24-34, 1983.) and WIGGLE (Lee, Sang-Ho, K. Palmo, S. Krimm (2005). WIGGLE: A new constrained molecular dynamics algorithm in Cartesian coordinates, Journal of Computational Physics 210: 171-182 2005). Both techniques are combined with the Verlet method or with one of the mentioned variants. Nevertheless, the success thereof in terms of increasing the interval of integration has been rather limited and have only provided factors of acceleration of between 2 and 5 (see Michael A. Sherman et al., Method for large timesteps in molecular modeling, U.S. patent Ser. No. 10/053,253), a completely insufficient improvement given the very long computing times provided, as mentioned, by the state of the art in molecular dynamics.
The most well known technique based on increasing the interval of integration of slow forces is known as MTS (multiple-timestep). In them, the conventional methods of integration in molecular dynamics (Verlet and its variants) are modified by evaluating the slow forces less frequently than the fast forces. The most popular and by far the most widespread MTS method in commercial simulators (GROMACS: H. J. C. Berendsen, D. van der Spoel and R. van Drunen, GROMACS: A message-passing parallel molecular dynamics implementation, Computer Physics Communications Volume 91, Issues 1-3, 2 Sep. 1995, Pages 43-56, AMBER: D. A. Pearlman, D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham, S. DeBolt, D. Ferguson, G. Seibel, P. Kollman, AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules, Computer Physics Communications Volume 91, Issues 1-3, 2 Sep. 1995, Pages 1-41, CHARMM: B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, M. Karplus, CHARMM: A program for macromolecular energy, minimization, and dynamics calculations, Journal of Computational Chemistry Volume 4, Issue 2, pages 187-217, Summer 1983, NAMD: J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chypot, R. D. Skeel, L. Kale, K. Schulten, Scalable molecular dynamics with NAMD, Journal of Computational Chemistry, 26:1781-1802, 2005,) is that known as r-Respa (M. E. Tuckerman, B. J. Berne, and G. J. Martyna, Reversible multiple time scale molecular dynamics. J. Chem. Phys., 97:1990-2001, 1992.), also called Verlet-I (see H. Grubmüller, H. Heller, A. Windemuth, and K. Schulten, Generalized Verlet algorithm for efficient molecular dynamics simulations with long-range interactions, Mol. Sim., 6:121-142, 1991) and which was proposed independently in the two preceding references. The philosophy of the r-Respa method is common for MTS methods, i.e., the less frequent evaluation of the fast forces, such that in each of these evaluations, the magnitude of the force which has not been taken into account between said evaluations is accumulated, these forces are thus updated by impulses. The improvement in the step of r-Respa integration is again marginal with typical times ranging from 0.5 fs, 2 fs, 4 fs (see H. Grubmüller, P. Tavan, Multiple Timestep Algorithms for Molecular Dynamics Simulations or f Proteins: How Good Are They?, Journal of Computational Chemistry, Vol. 19, No. 13, 1998) to 1 fs, 2 fs and 5 fs (see T. C. Bishop, R. D. Skeel, K. Schulten, Difficulties with Multiple Timestepping and Fast Multipole Algorithm in Molecular Dynamics, Journal of Computational Chemistry, Vol. 18, No. 14, 1997) and therefore very far from what is really needed in molecular dynamics. The main reason for this marginal r-Respa improvement is that the update of the slow forces by means of impulses happens periodically, particularly at a frequency close to multiples of some of the natural frequencies of the system. This circumstance means that the system enters into resonance (Schlick, M. Mandziuk, R. D. Skeel, and K. Srinivas, Nonlinear resonance artifacts in molecular dynamics simulations, J. Comput. Phys., 139:1-29, 1998) and the simulation is therefore not viable.
To prevent these resonances the family of integrators, MOLLY, (see B. Garcia-Archilla, J. M. Sanz-Serna & R. D. Skeel, Long-time-step methods for oscillatory differential equations, SIAM J. Sci. Comput. 20(1998), 930-963 and also J. M. Sanz-Serna, Mollified impulse methods for highly-oscillatory differential equations, SIAM J. Numer. Anal. 46(2008), 1040-1059) was introduced. Despite being an alternative to the use of restrictions, the timestep continues to be restricted to a few femtoseconds, allowing only a 50% increase in the timestep (J. A. Izaguirre, Q. Ma, T. Matthey, J. Willcock, T. Slabach, B. Moore, G. Viamontes, Overcoming instabilities in Verlet-I/r-RESPA with the mollified impulse method).
An alternative is to increase the mass of the hydrogen atoms because, due to their low mass, their bonds with other atoms are responsible for the higher oscillation frequencies. These methods allow increasing the timestep for the fast forces up to 4 fs without altering the dynamics, in comparison with the 2 fs of MOLLY, and reach up to 7 fs at the expense of eliminating the fastest oscillations of the system, but without modifying the thermodynamic equilibrium properties of the system (see K. A. Feenstra, B. Hess, H. J. C. Berendsen, Improving efficiency of large time-scale molecular dynamics simulations of hydrogen-rich Systems, Journal of Comput. Chem. Volume 20, Issue 8, pages 786-798, June 1999)
The molecular dynamics model will be seen below.
The basic unit of this model is the atom. A potential modeling the interactions between atoms is introduced for a system of N atoms.
The reason for introducing a potential instead of directly introducing the forces is both physical, because the energy can be written as the kinetic energy plus said potential, and mathematical, because potential is a scalar function.
The force can be easily obtained from the potential by simply applying the gradient, therefore if {right arrow over (r)}1, {right arrow over (r)}2, . . . , {right arrow over (r)}i, . . . , {right arrow over (r)}N are the coordinates of the N atoms, Newton's equation is written based on the potential UTOTAL as
                                          m            i                    ⁢                                                    d                2                            ⁢                                                r                  →                                i                                                    d              ⁢                                                          ⁢                              t                2                                                    =                              -                                          ∂                                                                                              ∂                                                      r                    →                                    i                                                              ⁢                                    U              TOTAL                        ⁡                          (                                                                    r                    →                                    1                                ,                                                      r                    →                                    2                                ,                …                ⁢                                                                  ,                                                      r                    →                                    N                                            )                                                          (        2        )            
Given the different types of occurring interactions, this potential is the sum of different contributions:UTOTAL=UBOND+UANGLE+UDIHEDRAL+UVDW+UCOULOMB  (3)
The more common analytical expressions thereof being the following
                              U          BOND                =                              ∑            bonds                    ⁢                                                    k                i                bond                            ⁡                              (                                                      r                    i                                    -                                      r                    i0                                                  )                                      2                                              (        4        )                                          U          ANGLE                =                              ∑            angles                    ⁢                                                    k                i                angle                            ⁡                              (                                                      θ                    i                                    -                                      θ                                          i                      ⁢                                                                                          ⁢                      0                                                                      )                                      2                                              (        5        )                                          U          DIHEDRAL                =                              ∑            dihedrals                    ⁢                                    ∑                              n                i                                      ⁢                                          k                i                dihed                            ⁡                              [                                  1                  +                                      cos                    ⁡                                          (                                                                                                    n                            i                                                    ⁢                                                      ϕ                            i                                                                          -                                                  γ                          i                                                                    )                                                                      ]                                                                        (        6        )                                          U          VDW                =                              ∑            i                    ⁢                                    ∑                              j                >                i                                      ⁢                          4              ⁢                                                ɛ                                      ij                    i                                                  ⁡                                  [                                                                                    (                                                                              σ                            ij                                                                                r                            ij                                                                          )                                            12                                        -                                                                  (                                                                              σ                            ij                                                                                r                            ij                                                                          )                                            6                                                        ]                                                                                        (        7        )                                          U          COULOMB                =                              ∑            i                    ⁢                                    ∑                              j                >                i                                      ⁢                                                            q                  j                                ⁢                                  q                  i                                                            4                ⁢                                  πɛ                  0                                ⁢                                  r                  ij                                                                                        (        8        )            
UBOND is called bond distance potential, ri being the distance at which the atoms forming the ith bond are located, ri0 is the equilibrium distance for said bond, and kibond is a constant quantifying the force of said bond, particularly of the harmonic potential used. This potential and the force derived therefrom model the fact that bound atoms are at a certain distance with respect to one another. For example, the distance of two carbon atoms bound by a single bond is about 1.54 amstrongs, which would be the equilibrium distance for this molecular model
UANGLE is called bond angle potential, θi being formed by 3 bound atoms, θi0 is the equilibrium angle for said angle, and kiangle is a constant quantifying the magnitude of said force, and which is related to the force and characteristics of the binding forces between the 3 atoms forming the angle. It models the angle formed by bound atoms. For example, for a water molecule, the equilibrium angle for the angle formed by the oxygen atom and hydrogen atom is about 105°.
UDIHEDRAL is called torsion or dihedral angle potential, φi being the angle formed by the planes defined by four bound atoms, γi is the equilibrium angle for said angle of torsion, ni is the so-called multiplicity, since given the periodic nature of this force various equilibrium angles can exist and kidihed is a constant quantifying said force. This potential models the planarity of certain molecules.
UVDW is the potential modeling the Van Der Waals type interactions between unbound atoms. The origin thereof is in the electromagnetic interactions between atoms due to the atomic structure, in this case the functional dependency is the so-called Lennard-Jones potential, where rij is the distance between the ith and jth atoms, and εij and σij are two constants quantifying the intensity of the interaction between that pair of atoms, having a long-range attraction nature and a very short-range repulsion nature.
UCOULOMB is the Coulomb potential modeling the electrostatic interactions between unbound atoms, where rij is the distance between the ith and jth atoms, qi and qj are their respective effective charges, and ε0 is the constant quantifying the intensity of the electrostatic force, called in this case the electric permittivity of the medium.
The various constants are obtained by semi-empirical means called force field constants (see J. W Ponder, D. A. Case, Force fields for protein simulations, Adv. Prot. Chem. 66 27-85 (2003)). They reflect the properties experimentally observed for the molecules. Given that the atomic bond and the atomic structure are governed by quantum mechanics, calculations based on quantum dynamics and on the experimental data are used for obtaining this Newtonian model of molecular dynamics.
Therefore, each force field tabulates the set of values necessary for modeling the molecule.
There are other possibilities for the functional dependency of the different potentials (see Chapter 8, Tamar Schlick, Molecular Modeling and Simulation: An Interdisciplinary Guide, Springer, 2006).
UBOND, UANGLE, UDIHEDRAL model the molecule structure due to the atomic bond. The effect thereof is to provide the experimental value for the distances and the angles between bound atoms. Given the force of the covalent bond, they are fast forces with high frequency oscillations.
The fact that they are calculated as a summation with respect to the number of bonds, angles and dihedral angles (torsions), the number of which is approximately proportional to N, is what determines the O(N) dependency of the computing cost.
UVDW and UCOULOMB in turn model the electromagnetic interactions between unbound atoms. The double summation is caused by the O(N2) dependency of the computing cost.
The effect of the periodic boundary conditions calculated by means of the aforementioned Particle Mesh Ewald method, providing an extra potential UEWALD, must be added to these potentials
For correct numerical integration of Newton's equation linking the coordinate variation with force and velocity, the rate of force variation is crucial.
This is because the greater the force variation between two time instants, the greater the numerical error.
A numerical integrator allows finding the coordinate value of the atoms in one instant based on a preceding instant, the time gap between both instants being the timestep h.
If the vector formed by the positions and velocities of the atoms of the system is called {right arrow over (q)}i, the following will be true:{right arrow over (q)}i+1={right arrow over (I)}(qi,h)  (9){right arrow over (I)}({right arrow over (q)}i,h) being the integrator. The coordinate dependency is essentially due to the force calculation, i.e.,{right arrow over (q)}i+1={right arrow over (I)}({right arrow over (F)}({right arrow over (q)}i),h)  (10)
Each integrator is characterized by a computing cost due to the force calculation FC and to the new coordinate calculation CC.
In the case of explicit integrators, this latter is significant. This is not the same for implicit integrators where this cost can be higher.
For good integration it is necessary that the force variation in interval h is small.
For the case of a variable timestep, predicting this variation is key in choosing the timestep.
For simulating the experimental conditions, the so-called thermostats (constant temperature) and barostats (constant pressure) are introduced, these conditions being present in actual systems. Berendsen's method (Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. (1984), Molecular-Dynamics with Coupling to an External Bath, Journal of Chemical Physics 81 (8): 3684-3690 (1984)), Nose-Hoover's method (see Nose, S, A unified formulation of the constant temperature molecular-dynamics methods, Journal of Chemical Physics 81 (1): 511-519 (1984) and also G J Martyna, D J Tobias, M L Klein, Constant pressure molecular dynamics algorithms, J. Chem. Phys 101(5), 1994) are commonly used, and the methods based on Langevin dynamics (see T. Schlick, Molecular Modeling and Simulation, Springer. pp. 435-438 (2002) and also S E Feller, Y. Zhang, R W Pastor and B R Brooks, Constant pressure molecular dynamics simulation: The Langevin piston method, J. Chem. Phys. 103(11), 1995) are particularly used. An integrator integrating Langevin dynamics is called a BBK integrator (see A. Brunger, C. L. Brooks, M. Karplus, Stochastic boundary conditions for molecular dynamics simulations of ST2 water, Chem. Phys. Lett., 105 (1984), pp. 495-500.)
The mentioned Langevin dynamics adds a stochastic term to Newton's equation, so energy is no longer conserved and the dynamics are not reversible, but it simulates the thermodynamic properties of a system at a constant temperature. A damping term proportional to velocity is further introduced. This term stabilizes the dynamics and allows increasing the timestep. However, in order for the dynamics to be realistic, the damping factor must not exceed a value, and the timestep reached must not exceed 14 fs for slow forces (see J. A. Izaguirre, D. P. Catarello, J. M. Wozniak, R. D. Skeel, Langevin stabilization of molecular dynamics, Journal of Chem. Phys. Volume 114 number 5 1 (2001)) or even 16 fs (see Q. Ma and J. A. Izaguirre, Targeted mollified impulse—a multiscale stochastic integrator for long molecular dynamics simulations). However, these methods use a high damping term which slows down the simulation, and therefore involves more iterations to reach the desired state, in addition to not reaching the timestep which would allow the dynamics of slow forces.
Despite recent computational advances which allow the millisecond simulation of a protein in a reasonable time (see D. E. Shaw, R. O. Dror, J. K. Salmon, J. P. Grossman, K. M. Mackenzie, J. A. Bank, C. Young, M. M. D., B. Batson, K. J. B., E. Chow, M. P. Eastwood, D. J. Ierardi, J. L. Klepeis, J. S. Kuskin, R. H. Larson, K. Lindorff-Larsen, P. Maragakis, M. A. Moraes, S. Piana, Y. Shan, B. Towles, Millisecond-Scale Molecular Dynamics Simulations on Anton, Proceedings of the Conference on High Performance Computing, Networking, Storage and Analysis (SC09), New York, N.Y.: ACM Press, 2009), these are based on solely computational improvements. They do not however overcome the femtosecond scale for the integration step, so there is still a need for enormous computational resources that are beyond the reach of many.